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Abstract 

This  thesis  addresses  the  problem  of  determining  aircraft  position  during 
flight  given  noisy  and  biased  measurements  from  a  barometric  altimeter  and  two 
tactical  air  navigation  (TACAN)  transceivers. 

A  Kalman  smoother  is  developed  to  perform  post-flight  data  processing  on 
the  measurement  data.  The  smoother  estimates  aircraft  position,  velocity,  and 
acceleration  as  well  as  biases  in  the  measurements. 

Since  actual  flight  test  data  is  not  available,  computer  simulations  examine 
the  performance  of  this  tracking  technique.  The  simulated  flight  includes  low  and 
high-speed  turns,  constant  rate  ascents,  descents,  and  accelerations.  The  tracking 
algorithm  tracked  the  aircraft  to  within  135  meters  of  its  actual  position  95  percent 
of  the  time.  An  unaided  inertial  navigation  system  used  by  the  4950th  Test  Wing 
in  another  flight  test  program  showed  a  position  error  growth  rate  of  2000  meters 
per  hour. 

The  computer  programs  which  perform  the  smoothing  are  written  in  Fortran- 
77  and  run  under  the  Digital  Equipment  Corporation  VMS  operating  system. 


AIRCRAFT  TRACKING  WITH  DUAL  TACAN 


I.  Introduction 


1.1  Overview 

The  4950th  Teat  Wing,  located  at  Wright-Patteraon  AFB,  teats  new  avionics 
and  electronic  warfare  equipment  for  the  Air  Force.  Performance  of  these  systems 
often  varies  as  the  aircraft  position  changes.  While  recording  test  item  data  during 
a  flight  test,  the  Wing  also  records  Time/Space  Position  Information  (TSPI),  a 
record  of  the  airplane’s  flight  path.  TSPI  allows  Wing  data  analysts  to  account  for 
a  system’s  varying  performance  under  test. 

The  airborne  navigation  and  tracking  systems  currently  available  in  the 
4950th  Test  Wing  do  not  provide  sufficiently  accurate  TSPI  for  evaluating  test 
item  performance.  This  inability  to  record  flight  path  information  forces  the  Test 
Wing  to  fly  test  missions  at  ranges  instrumented  with  radars  to  track  the  airplane. 
Since  Wright- Patterson  AFB  does  not  have  an  instrumented  test  range,  the  Wing 
deploys  aircraft  and  personnel  to  other  bases  for  test  flights. 

Deploying  aircraft  and  personnel  to  test  ranges  has  many  drawbacks,  includ¬ 
ing  costs  for  travel,  logistical  support,  and  range  operating.  These  additions  can 
easily  triple  the  cost  of  collecting  flight  test  data.  Deploying  to  a  test  range  can 
also  adversely  affect  the  program  schedule.  Conflicts  over  use  of  range  resources 
arise  when  programs  from  different  test  organizations  both  require  use  of  range 
facilities.  If  the  4950th  Test  Wing’s  program  has  a  low  Department  of  Defense 
priority,  it  may  have  to  wait  for  an  opportunity  to  fly  on  the  range. 
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l.t  Problem  Statement 


The  4950th  Test  Wing  needs  an  inexpensive  and  accurate  aircraft  tracking 
system.  The  Wing  wants  to  know  if  optimally  combining  the  signals  from  two  Tac¬ 
tical  Air  Navigation  (TACAN)  transceivers  can  supply  this  tracking  information 
[21].  This  will  provide  the  capability  to  fly  test  missions  without  using  a  test  range 
radar  for  tracking  support. 

1.3  Assumptions 

The  Test  Wing  does  not  need  tracking  data  displayed  in  real  time  on  the 
aircraft.  Post  flight  data  analysis  allows  the  use  of  existing  Test  Wing  computer 
facilities.  A  real-time  system  requires  installation  of  a  computer  on  the  aircraft, 
increasing  system  cost  and  complexity. 

The  aircrew  selects  TACAN  beacons  to  provide  coverage  during  each  flight 
segment.  They  decide  which  TACAN  beacons  to  use  while  mission  planning  before 
the  test  flight. 

The  airborne  system  records  station  identifiers,  ranges,  and  bearings  from 
TACAN  beacons.  The  system  also  records  the  time  of  day  and  barometric  altitude. 
This  information  allows  post-flight  data  processing  to  calculate  a  complete  three- 
dimensional  time  history  of  the  aircraft  flight  path. 

1.4  Standards 

The  4950th  Test  Wing  requires  TSPI  accurate  to  within  several  hundred  feet 
of  the  actual  aircraft  position.  In  addition  to  meeting  this  accuracy  requirement., 
the  tracking  system  employed  must  conduct  any  necessary  pre-flight  calibrations 
rapidly.  The  delay  between  equipment  power-up  and  takeoff  must  be  less  than  30 
minutes.  Converting  TACAN  measurement  files  into  aircraft  track  files  should  not 
slow  the  analysis  of  project  data.  Less  than  two  days  of  data  reduction  should 
provide  TSPI  data  for  a  two-hour  flight. 
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The  track  files  generated  by  this  system  can  be  formatted  to  match  TSP1 
files  generated  at  Eglin  AFB.  This  compatibility  avoids  the  need  to  develop  two 
sets  of  analysis  tools;  one  for  use  with  dual  TACAN  data  and  another  for  use  with 
test  range  radar  data.  Test  programs  would  thus  have  the  flexibility  of  flying  with 
either  source  of  tracking  data. 

Since  the  dual  TACAN  system  cannot  provide  extremely  accurate  tracking 
data,  it  would  not  replace  all  test  range  flying.  Dual  TACAN  can  provide  a  low- 
cost  alternative  for  programs  not  requiring  extremely  accurate  positioning  data. 
Dual  TACAN  would  also  allow  local  flying  to  verify  test  system  operation  before 
deploying  a  test  team  to  an  instrumented  range  for  more  comprehensive  flight  tests. 

If  the  Dual  TACAN  system  cannot  be  readily  implemented,  the  Wing  will 
continue  to  use  their  current  limited  tracking  ability  until  the  Global  Positioning 
System  (GPS)  becomes  operational  and  user  equipment  becomes  available.  When 
fully  deployed,  GPS  satellites  will  allow  an  aircraft  to  record  its  position  within  15 
meters  [19:55].  Dual  TACAN  can  serve  as  an  interim  measure  until  GPS  becomes 
fully  operational  and  the  Test  Wing  installs  GPS  receivers  on  their  aircraft. 

1.5  Scope 

This  project  builds  upon  work  already  performed  in  the  Test  Analysis  Branch 
of  the  4950th  Test  Wing.  They  have  developed  a  program  to  take  TACAN  mea¬ 
surement  data  and  convert  it  into  an  aircraft  track  file.  Because  of  budget  and 
manpower  limitations,  their  programs  do  not  attempt  to  filter  out  noise  from  data. 
Digitally  processing  TACAN  information  with  a  Kalman  filter  can  remove  noises 
and  biases  to  produce  a  more  accurate  track  file. 

The  proposed  system  uses  two  AN/ARN-118(V)  TACAN  transceivers.  The 
AN/ARN-118  serves  as  the  standard  TACAN  system  for  Test  Wing  aircraft  [39:3]. 
Taking  additional  measurements  from  a  third  or  fourth  TACAN  beacon  by  in¬ 
stalling  additional  transceivers  or  continuously  retuning  the  available  transceivers 
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can  lead  to  a  more  accurate  tracking  system.  These  additional  measurements  are 
not  considered  in  this  project  because: 

1.  Installing  additional  TACAN  transceivers  is  difficult  because  of: 

•  limited  available  airplane  equipment  space 

•  limited  available  airplane  equipment  power 

•  requirement  to  mount  additional  TACAN  antennas  externally 

2.  Manually  retuning  available  transceivers  increases  the  aircrew  workload  and 
has  the  potential  for  errors  and  loss  of  track  data. 

3.  Automatically  retuning  TACAN  transceivers  represent  costly  custom  equip¬ 
ment  not  currently  available  in  the  Test  Wing. 

4.  Constantly  retuning  the  airborne  transceivers  forces  the  TACAN  to  operate 
in  an  acquisition  mode  which  may  last  as  long  as  20  seconds.  During  this 
time,  the  transceiver  places  an  operating  load  6  times  the  standard  tracking 
level  on  the  ground  beacon.  This  can  seriously  degrade  the  beacons  ability 
to  handle  multiple  aircraft  [7:15]. 

1.6  Approach 

The  initial  task  in  this  project  involves  reviewing  technical  literature  on  nav¬ 
igation  and  signal  processing  systems.  This  provided  background  on  how  a  dual 
TACAN  system  can  provide  tracking  data  and  benefits  possible  from  filtering  that 
data.  Chapter  II  of  this  thesis  summarizes  the  results  of  the  literature  review. 

In  1986  and  1987,  the  4950th  Test  Wing  performed  preliminary  work  with 
TACAN  tracking.  While  a  lack  of  funding  and  manpower  prevented  completion  of 
the  project,  the  Wing’s  Test  Analysis  Division  (4950  TESTW/FFT)  used  computer 
simulations  to  analyze  a  TACAN  tracking  system’s  performance  [39,17],  These 
programs  are  the  starting  point  for  this  thesis. 


This  project  does  not  involve  any  special  equipment.  AFIT’s  computer 
systems  provide  the  development  environment  for  new  programs  which  convert 
TACAN  measurement  files  into  TSPI  track  files.  Program  checkout  and  test  data 
analysis  use  both  AFIT  and  Test  Wing  computers. 

This  analysis  requires  an  accurate  model  for  the  TACAN  transceiver  mea¬ 
surements.  Chapter  III  develops  a  model  for  the  tracking  problem.  Chapter  IV 
concentrates  on  the  actual  design  and  development  of  the  tracking  software.  Chap¬ 
ter  V  discusses  the  simulations  conducted  using  the  dual  TACAN  tracker  and  com¬ 
pares  the  accuracy  of  the  filter  reported  position  against  actual  aircraft  position. 
Chapter  VI  summarizes  the  project. 
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II.  Background 


2.1  Overview 

Flight  testing  new  avionics  and  electronic  warfare  equipment  for  the  Air 
Force  usually  requires  accurate  knowledge  of  the  aircraft’s  position  throughout  the 
flight.  Currently  available  air  navigation  systems  in  the  4950th  Test  Wing  cannot 
provide  this  type  of  information.  This  forces  the  Wing  to  fly  their  test  missions  at 
instrumented  test  ranges,  such  as  Eglin  AFB,  FL  or  Holloman  AFB,  NM. 

The  Test  Wing  would  like  to  have  a  system  on  the  aircraft  capable  of  pro¬ 
viding  this  positioning  data.  This  system  would  allow  projects  to  conduct  tests 
without  deploying  to  test  ranges. 

2.2  Navigation  and  Tracking  Systems 

The  Test  Wing  currently  has  several  different  navigation  systems  available  to 
support  test  programs.  All  Wing  aircraft  are  equipped  with  TACAN  transceivers. 
Standard  TACAN  measurements  are  not  very  accurate,  and  position  uncertainties 
increase  with  aircraft  range  from  the  ground  beacon. 

The  Test  Wing  also  uses  Carousel  IV-E  inertial  navigation  systems  (INS). 
These  are  not  high-precision  INS  units  by  current  standards.  The  Carousel  IV- 
E  accuracy  specification  allows  the  radial  position  error  to  grow  at  a  rate  of  two 
nautical  miles  per  hour  [39:26]. 

The  satellite- based  GPS  is  the  newest  system  available  in  the  Test  Wing.  In 
the  near  future,  the  flexibility  and  high  accuracy  of  GPS  will  satisfy  virtually  all  of 
the  Wing’s  aircraft  tracking  requirements.  Unfortunately,  airborne  GPS  equipment 
is  not  yet  readily  available,  and  all  the  satellites  necessary  to  support  GPS  have 
not  been  launched  [19]. 


As  an  interim  measure,  John  Franzen,  a  data  analyst  in  the  4950th  Test 
Wing,  investigated  methods  to  process  available  navigation  data  after  a  flight  to 
obtain  more  accurate  estimates  of  the  aircraft  track  [15,17,18].  A  Test  Wing  report 
suggested  a  TACAN-based  system  to  provide  flight  path  data  [39]. 

2.3  TA  CAN  Principles 

Small  size,  low  power  requirements,  simple  operation,  and  low  cost  all  com¬ 
bine  to  make  TACAN  the  primary  method  of  aircraft  navigation  in  the  United 
States  [19:50].  While  the  Department  of  Defense  does  operate  some  TACAN  sites, 
the  majority  of  beacons  installed  in  the  United  States  are  maintained  by  the  Federal 
Aviation  Administration,  a  branch  of  the  Department  of  Transportation. 

The  airborne  TACAN  equipment  provides  two  pieces  of  information  to  the 
crew.  By  transmitting  pulses  to  a  ground  station  and  timing  the  delay  before 
receiving  response  pulses,  the  system  measures  distance  between  the  aircraft  and 
beacon.  The  system  also  displays  bearing  from  the  beacon  to  the  aircraft  based  on 
variations  in  the  beacon’s  rotating  antenna  pattern. 

Navigation  with  TACAN  involves  polar  geometry.  The  airborne  equipment 
displays  the  aircraft  range  and  bearing  from  a  ground  beacon.  Beacon  locations 
are  displayed  on  navigation  charts,  making  it  simple  to  determine  the  airplane’s 
position  with  a  chart,  a  pair  of  dividers  and  a  protractor  (see  Figure  1). 

To  simplify  enroute  navigation  using  TACAN,  the  Federal  Aviation  Admin- 
istraton  has  established  airways  and  jet  routes  between  TACAN  beacons.  Aircraft 
fly  along  these  routes  from  one  TACAN  to  another,  using  the  range  measurement 
to  indicate  their  position  in  the  route  and  providing  an  estimate  of  how  long  until 
the  next  course  change. 
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g-4  TACAN  Shortcomings 

Polar  geometry  causes  TACAN  accuracy  to  decrease  as  distance  to  the  station 
increases.  The  magnitude  of  typical  TACAN  measurement  errors  average  1.65 
degrees  of  bearing  error  and  600  feet  of  range  error  [39:50].  For  an  aircraft  located 
30  miles  from  a  TACAN  beacon,  1.65  degrees  of  bearing  error  and  600  feet  of  range 
error  represent  over  5000  feet  of  position  error  (see  Figure  2). 

Inaccurate  bearing  information  comprises  the  major  error  component  of 
TACAN  signals.  Combining  several  more  accurate  range  measurements  can  de¬ 
termine  an  aircraft’s  position  without  using  bearing  data.  Groginsky  showed  how 
range  measurements  to  three  different  ground  beacons  could  provide  the  position 
and  altitude  of  an  aircraft  [22:178].  Figure  3  shows  the  ideal  case  of  three  range 
measurements  providing  a  position  fix.  Measurement  noises  and  biases  cause  the 
typical  situation  to  resemble  Figure  4  with  no  single  unambiguous  aircraft  position 
indicated  by  the  measurements. 
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Figure  2.  Single  TACAN  Positioning  Error 


Figure  3.  Ideal  Triple  Range  Fix 


2.5  Multiple  TA  CAN  Processing 

Multiple  range  measurements  for  aircraft  positioning  served  as  the  basis  for 
Latham’s  research  at  Grumman  Aerospace  Corporation.  In  a  prototype  system, 

l 

Latham  combined  range  data  from  ten  different  ground  beacons  and  could  consis¬ 
tently  calculate  an  aircraft’s  position  within  200  feet  of  the  actual  position  [26:157]. 

Later  work  added  inertial  navigation  system  data  for  greater  accuracy  [27:72]. 

Riggins  followed  on  Latham’s  work  and  examined  different  methods  of  com¬ 
bining  multiple  range  measurements  [47].  The  technique  developed  by  Riggins 
could  calculate  aircraft  position  during  flight  without  need  for  post-flight  com¬ 
puter  analysis  of  the  range  data.  This  system  accurately  positioned  the  aircraft  * 

* 

within  200  feet  [47:101]. 

Latham  and  Riggins  used  single-channel  TACAN  systems  on-board  the  air¬ 
craft.  To  obtain  multiple  range  measurements,  computers  automatically  switched  ^ 

the  aircraft  TACAN  transceiver  to  different  ground  beacons  throughout  the  flight. 

Existing  Air  Force  TACAN  equipment  used  in  the  Test  Wing’s  aircraft  cannot 

perform  automatic  switching.  { 
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I 

Tests  conducted  by  the  4950th  Test  Wing  suggest  a  different  approach  to 
i  position  fixing.  The  Wing’s  proposed  system  would  record  the  output  from 

two  separate  TACAN  transceivers  along  with  the  output  of  a  pressure  altime¬ 
ter.  These  three  measurements  provide  sufficient  information  to  determine  aircraft 
|  position  [44:28]. 

Analysis  of  TACAN  signals  by  the  Test  Wing  revealed  the  potential  to  ex¬ 

tract,  with  additional  signal  processing,  high  accuracy  range  measurements  by 
j  removing  measurement  biases.  Error  analysis  based  on  measurements  of  TACAN 

stations  near  Wright- Patterson  AFB  indicate  that  such  a  system  could  accurately 
position  the  airplane  within  300  feet  [39:F7].  After  initial  analysis  of  the  problem, 
lack  of  manpower  and  funding  prevented  the  Wing  from  investigating  this  system 
^  further. 

Systems  currently  in  development  offer  tracking  capability  far  exceeding  the 
Test  Wing’s  requirements.  The  satellite-based  GPS  will  provide  position  infor¬ 
mation  accurate  to  within  15  meters  (19:55j.  Although  GPS  was  scheduled  to  be 
fully  operational  by  1984  [19:55],  the  seven  satellites  currently  in  orbit  are  pre- 
production  models.  The  first  production  satellite  will  be  launched  30  December 
1988  with  6  additional  satellites  added  each  year  for  the  next  three  years  [20].  At 
this  rate,  the  21  satellites  necessary  for  full  GPS  coverage  will  not  be  available  un¬ 
til  1991.  Even  if  the  satellites  were  orbiting,  airborne  GPS  receivers  are  currently 
expensive  and  difficult  to  obtain. 

2.6  Summary 

Combining  the  signals  from  multiple  TACAN  transceivers  offers  an  oppor¬ 
tunity  to  collect  precision  tracking  data  without  the  need  for  special  equipment. 
This  will  improve  the  Test  Wing’s  in-house  tracking  capability  in  an  economical 
manner. 
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III.  Model  Development 


The  system  model  serves  as  the  heart  of  any  stochastic  estimator.  The  model 
outlined  in  this  chapter  is  driven  by  the  available  measuring  devices  and  a  simplified 
model  of  the  aircraft  dynamics. 

3.1  Choice  of  States 

On-board  navigation  systems  typically  access  commands  entered  by  the  pilot 
and  use  this  information  when  estimating  aircraft  trajectory.  The  T-39B  Sabre- 
liners  flown  by  the  4950th  Test  Wing  have  no  provisions  for  measuring  or  record¬ 
ing  this  information.  Without  command  inputs,  the  tracking  system  resembles  a 
weapon  guidance  system  tracking  an  uncooperative  target. 

Aircraft  dynamics  were  divided  into  three  orthogonal,  uncoupled  Cartesian 
components.  To  represent  unknown  pilot  inputs,  the  model  uses  three  independent 
white  Gaussian  noises  arranged  in  a  vector  denoted  as  w (t).  To  represent  physi¬ 
cal  maneuvering  characteristics  of  the  aircraft  more  accurately,  this  input  passes 
through  a  first-order  lag  shaping  filter  with  a  time  constant  of  T  to  correlate  input 
accelerations.  Integrating  acceleration  provides  velocity;  integrating  velocity  yields 
position.  This  type  of  dynamics  model  is  often  used  for  non-cooperative  aircraft 
in  targeting  systems  [6,38,47]. 

Equations  (l)-(3)  state  the  scalar  version  of  this  model  for  one  dimension: 


P(t) 

=  V(t) 

(1) 

V(t) 

=  A(t) 

(2) 

A(t) 

=  ~~A{t)  -1-  w{t) 

(3) 

Repeating  the  equations  in  the  x ,  y  and  z  directions  provide  nine  separate 
equations,  driven  by  three  independent  noise  sources. 
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Figure  5.  Earth- Centered,  Eartli-Fixed  Coordinate  System 


The  orthogonal  xyz  reference  frame  uses  earth-centered,  earth-fixed  coordi¬ 
nates.  This  places  the  origin  at  the  earth’s  center  with  the  i-axis  pointing  out 
along  the  equator  at  the  prime  meridian,  the  z-axis  pointing  through  the  north 
pole,  and  the  t/-axis  pointing  out  along  the  equator  at  90  degrees  east  longitude  to 
complete  a  right-handed  xyz  coordinate  system.  Figure  5  shows  the  earth-centered, 
earth-fixed  coordinate  system. 

Expressing  the  position,  velocity  and  acceleration  components  as  a  vector 
differential  equation  models  the  aircraft  dynamics  equations  in  standard  form  as  a 


ninth-order  system  (with  r  representing  —1/T  to  allow  more  compact  equations): 
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(4) 


The  value  — (l/4)sec  1)  used  for  r  corresponds  to  the  acceleration  time  con¬ 


stant  found  in  a  typical  non-fighter  aircraft  [37]. 


For  navigation  purposes,  the  earth  is  modeled  as  an  ellipsoid  with  semi-major 
axis  a  and  eccentricity  e  [1],  Based  upon  this  reference  ellipsoid,  transforming 
geodetic  coordinates  expressed  in  latitude,  longitude  and  altitude  (L,  1,  h)  to  earth- 
centered,  earth-fixed  coordinates  used  by  the  filter  ( Px ,  Pv,  Pz)  is  done  with: 


P r 
Py 
P r 


,  +  h 

Vl  —  €2  sin*  L 

_ a _ ^ 

Vl  —  e2  sin2  L 

Vl  -  eJ  sin2  L 


cos  L  cos  l 
cos  L  sin  l 
sin  L 


(5) 

(6) 
(7) 


In  addition  to  the  nine  states  already  selected,  the  filter  must  model  char¬ 
acteristics  of  the  measurement  sensors.  One  measurement  available  to  the  filter  is 
pressure  altitude.  Delay  in  the  altitude  transducer  requires  the  addition  of  another 
state  to  the  filter  which  accounts  for  the  lag  in  altimeter  output  during  ascents 
and  descents.  The  altimeter  value  depends  on  the  previous  altimeter  reading  and 
the  actual  altitude.  Altitude  is  a  nonlinear  function  of  the  aircraft  position  vector 
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P  consisting  of  the  components  Pat  Pv,  Pt.  This  function  gives  the  length  of  the 
vector  from  the  aircraft  position  perpendicular  to  the  reference  geoid. 


MO  =  -jK(  0  +  ^MP(0) 

=  aha(t )  —  ah{P(t)} 


(8) 


Equation  (9)  shows  the  nonlinear  relationship  between  the  aircraft  earth- 
centered,  earth-fixed  position  and  altitude. 


P.  a(l  -  e1) 
sin  L  y/\  —  «2  sin*  L 


(9) 


This  nonlinearity  forces  the  tracker  to  use  an  extended  Kalman  filter. 


Along  with  altimeter  dynamics,  all  measurements  contain  biases  the  filter 
must  estimate  and  remove.  Tests  conducted  by  the  4950th  Test  Wing  and  by 
the  British  Royal  Aircraft  Establishment  indicate  TACAN  errors  consist  of  large, 
slowly  varying  biases  ranging  from  —225  to  +340  meters  and  small  white  noise 
components  with  approximate  covariance  of  400m2  [39,45].  This  tracking  system 
uses  five  measurements,  each  with  an  unknown  bias.  These  result  in  filter  states 
for  the  two  range  biases  (RBX,  RB2),  the  two  bearing  biases  (BBX,  BB2)  and  the 
altimeter  bias  ( hB ).  Each  bias  state  corresponds  to  the  output  of  an  integrator 
driven  by  small  magnitude,  pair-wise  uncorrelated  noises.  The  initial  values  of 
these  integrators  are  read  from  a  data  file  based  on  previous  flight  tests. 
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(10) 


When  combined,  the  preceding  equations  describe  a  system  using  the  fifteen 
states  shown  in  Table  1.  Since  the  bias  states  are  independent  of  position,  veloc- 
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Table  1.  Tracking  System  States 


H 

Aircraft  Position 

t§ 

Aircraft  Velocity 

i 

Aircraft  Acceleration 

H2H 

Altimeter  Output 

RBi 

BBx 

TACAN  1  Biases 

rb2 

bb2 

TACAN  2  Biases 

mam 

Altimeter  Bias 

ity,  acceleration  or  altitude,  bias  states  can  be  appended  to  provide  the  following 
system  representation: 


Px 

— 

(11) 

A 

= 

Vy 

(12) 

A 

= 

K 

(13) 

Vm 

= 

Ax 

(14) 

Vy 

Ay 

(15) 

vz 

= 

Az 

(16) 

Ax 

= 

1  . 

-  j,Ax  +  wx 

(17) 

K 

- 

1  . 

-fAy+Wy 

(18) 

At 

= 

1  . 

~jAz  +  wz 

(19) 

K 

= 

_  jT^O  +  jTh{PX,  Py,  Pz) 

(20) 

RBt 

— 

WRBi 

(21) 
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Table  2.  Available  Measurements 


Ri 
B x 

slant  range  to  beacon  1 

magnetic  heading  from  beacon  1  to  aircraft 

R2 

b2 

slant  range  to  beacon  2 

magnetic  heading  from  beacon  2  to  aircraft 

h? 

pressure  altitude  measured  on  aircraft 

BB\  =  wbbi 

(22) 

RB2  —  WRBi 

(23) 

BB2  =  wbbi 

(24) 

hB  =  WhB 

(25) 

This  can  be  stated  more  succinctly  in  matrix  form  as: 


x(<)  =  f[x(()]  +  G(i)w(() 


(26) 


3.2  Measurement  Model 

Having  established  the  dynamics  model  for  the  target  aircraft,  the  next  step 
in  designing  a  filter  is  to  determine  the  measurement  model.  The  available  mea¬ 
surements  shown  in  Table  2  can  be  related  to  the  fifteen  system  state  variables  to 
establish  a  measurement  model  for  the  state  estimator. 

The  following  equations  show  the  relationship  between  state  variables,  re¬ 
corded  sampled-data  measurements,  and  the  components  of  the  noise  vector  v(tj) 
entering  each  measurement.  The  functions  for  range  and  bearing  are  discussed 
later  in  this  chapter. 


R\{Px  j  Py )  Pt)  +  RB i  +  l>Ri 

(27) 

B\{Px,Py,Pz )  +  BBi  4-  vb i 

(28) 

Rt{  Px )  Py,  Pz)  +  RB2  4  VR2 

(29) 
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(30) 

(31) 


i?j  =  5j(P*,  Py ,  PZ)  -f  BBl  +  VBt 
hp  =  ha  +  hB  +  vhp 

This  set  of  equations  can  be  rewritten  in  vector  form  as: 

z(f,)h(x(ti)}  +  y(U)  (32) 

The  filter  knows  a  priori  the  earth-centered,  earth- fixed  coordinates  for  the 
TACAN  beacons  (XB,  Yb  and  Zb)-  Using  TACAN  position  and  aircraft  coordi¬ 
nates,  the  following  equation  gives  slant  range  from  the  beacon  to  the  aircraft: 

R  =  \J ( Px  -  XB)2  +  (Pv  -YBy  4-  ( Pz  -  ZB)2  (33) 

The  filter  also  uses  bearing  measurements  from  the  TACANs  to  the  air¬ 
craft.  These  measurements  are  magnetic  bearings  based  on  the  local  tangent  plane 
defined  at  the  beacon.  The  range  vector  between  the  beacon  and  the  aircraft, 
[(P*  —  Xb)  (Pv  —  Yb)  (Pz  —  Zb)]T,  is  projected  into  a  north,  east,  down  coordi¬ 
nate  system  based  on  the  beacon’s  geodetic  latitude  (Lb)  and  longitude  (lB)-  The 
north,  east  and  down  (N,  E,  D)  components  of  the  range  vector  are  given  by  the 
vector  transformation  [4]: 

N  —  sin  Lb  cos  Ib  —sin  Lb  sin  Ib  cos  Lb  Px  —  Xb 

E  =  —  sin  Ib  cos  lB  0  Py  —  YB  (34) 

D  ~cosLBcoslB  —  cos  Lb  sin /b  —sin  Lb  Pz  —  ZB 

After  calculating  the  north  and  east  range  vector  components,  using  the  in¬ 
verse  tangent  gives  true  heading  from  the  beacon  as  shown  in  Figure  6.  The  mag¬ 
netic  heading  reported  by  the  TACAN  receiver  will  also  have  the  local  magnetic 
variation  at  the  beacon  added  to  the  true  heading.  Note  that  east  magnetic  vari¬ 
ation  is  added  while  west  magnetic  variation  is  subtracted  from  the  true  heading, 
(3,  to  give  magnetic  heading, B.) 
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(3  =  tan'1  |  (35) 

B  =  (3  +  var  (36) 


The  transformation  matrix  defined  in  Equation  (34)  depends  only  on  beacon 
coordinates,  not  aircraft  position.  A  lookup  table  stores  the  geographic  coordinates 
and  magnetic  variation  for  each  TACAN  beacon. 


3.3  Summary 

The  15-state  dynamics  model  and  the  5-state  measurement  model  discussed 
in  this  chapter  serve  as  the  foundation  for  developing  a  state  estimator.  The  next 
chapter  develops  a  Kalman  filter  to  process  the  measurement  data  and  estimate 
the  aircraft  position. 
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IV.  Design  and  Development 


The  Chapter  III  measurement  and  dynamics  models  serve  as  the  basis  for 
a  stochastic  estimator  which  takes  noisy  measurements  and  estimates  the  aircraft 
flight  path. 

4-1  Kalman  Filters 

The  tracking  system  developed  in  this  study  uses  a  Kalman  filter  to  estimate 
errors  in  sensor  measurements  and  to  allow  accurate  calculation  of  aircraft  position. 
A  Kalman  filter  is  an  optimal  recursive  data  processing  algorithm  designed  to 
estimate  state  variables  of  interest.  The  filter  uses  all  available  measurements, 
regardless  of  their  precision,  as  well  as  a  priori  knowledge  of  system  dynamics 
[34:4]. 

4.1.1  Linear  Systems  The  basic  equations  for  a  linear  Kalman  filter  involve 
the  state  dynamics,  the  solution  to  a  stochastic  difference  equation,  a  measurement 
model,  and  the  equations  for  a  measurement  update.  Equations  (37)— (46)  show 
these  relationships  in  standard  matrix  notation.  The  dynamics  model  for  the 
system  is: 

x(t)  =  F  (t)x(t)  +  B(t)u(t)  +  G(t)w(t)  (37) 

with  the  system  state  vector  x(t),  the  control  input  u(<)?  and  a  zero-mean  white 
Gaussian  noise  w(t).  This  noise  is  assumed  independent  of  x(t)  for  all  time,  and 
has  a  strength,  Q(t),  given  by: 

E{w(t)wT{t  +  r)]  =  Q(f)«(r)  (38) 

The  filter  propagation  cycle  between  measurements  is  given  by: 

*{*n  =  #(<„<i-1)x(C_x)+  T  #(«,-, r)B(r)u(r)dr  (39) 
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p  (<f)  =  *(<;,*i-l)p(*t1)*r(Mi-l) 

+  r  #(<i,r)G(r)Q(T)Gr(r)#r(<i,<i.1)<iT  (40) 

where  $(f,r)  is  the  state  transition  matrix  associated  with  F(t).  The  state  tran¬ 
sition  matrix  solves  the  differential  equation: 

*(tfT)  =  F(t)*(t,T)  (41) 

subject  to  the  constraint 

$(t,t)  =  I  (42) 

The  discrete  time  measurements  available  to  the  system,  z(/,),  are  corrupted 
by  zero-mean  white  Gaussian  discrete-time  noise  v(f,)  which  has  a  covariance  of 
R(tj).  The  measurement  noise  v(f)  is  assumed  independent  of  x(t)  and  w (t)  for 
all  time.  The  measurement  are  related  to  the  system  state  vector  by: 


z(<)  =  H(<i)x(fj)  +  v(U)  (43) 

The  measurement  update  cycle  for  the  filter  is: 

K(«<)  =  P(tr)HT(ti)  [H(f,)P(f-)Hr(ft)  +R(<,)]_1  (44) 

*(tf)  =  x(t7)  +  K(ti)[z(ti)-H(ti)x{t~)]  (45) 

P (K)  =  P(*f )  “  K(<i)H(<i)P(<r )  (46) 


In  the  Kalman  filter  equations,  x  represents  the  filter’s  best  estimate  of  the 
actual  system  state  x,  P  represents  the  state  error  covariance  matrix,  and  the 
vector  z  represents  measurements  available  to  the  filter.  The  time  argument  t~ 
represents  the  instant  before  incorporating  a  measurement,  tf  the  instant  after 
using  a  measurement  to  update  the  system.  For  a  constant  F  matrix,  the  state 
transition  matrix  $  is  given  by: 

=  £  1  {i,I-Fi  <47> 
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4-l.t  Nonlinear  Systems  For  this  aircraft  tracking  problem,  the  propaga¬ 
tion  and  measurement  equations  used  in  the  system  model  are  nonlinear  functions 
of  aircraft  position.  This  makes  it  necessary  to  use  an  extended  Kalman  filter.  The 
dynamics  model  for  this  formulation  is: 

x(t)  =  f[x(t),  u(t),  t\  +  G(f  )w(< )  (48) 

The  propagation  equations  for  this  system  are: 

*(<D  =  /’  f(x(r),u(r),T]dT  (49) 

Jti- 1 

P  (if)  =  ^[ft>  U-i ;  x(r  )]P(t»-i  )#[<<,  <i_i ;  x(r )] 

+  f  !<_!-,  x(r)]G(r)Q(T)Gr(r)$7’[fi,f<_1;x(r)]  dr  (50) 

Jti-i 

The  function  $[f,-,f,_i;x(r)]  represents  the  solution  to  the  differential  equa¬ 


tion: 

*[<i,«<_i,x(r)]  =  F[«;x(r)]*[<j,<i_i,x(r))  (51) 

* 

with  the  boundary  condition: 

*[ti,ti,x(r)]  *=  I  (52) 

The  measurement  model  and  update  equation  for  the  nonlinear  filter  are: 
z(ti)  =  h[x(ii),u(<i)]  +  v(f,-)  (53) 

K(ti)  =  P(fr)Hr[f<;x(f-)i{H(f,;x(fr)JP(t-)HT[<i;x(f-)]  +  R(f,)}"1  (54) 

x{tt)  =  x{tr)  +  K(fi){z(f<)  -  h(x(<, '),«<]}  (55) 

P  (if)  =  P(<r)-K(f,)H[fi;x(f-)]P(<r)  (56) 


The  extended  Kalman  filter  recalculates  the  elements  of  the  F  and  H  matrices 
during  each  update  and  propagation  cycle  by  evaluating  partial  derivatives  of  the 
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vector  functions  f  and  h: 


H[<i;x(<)] 


*=*(«." ) 

5h(x,t) 


(57) 

(58) 


When  applying  an  extended  Kalman  filter  to  a  real-time  application,  the 
requirement  to  recalculate  F  and  H  continually  can  impose  a  severe  computational 
burden.  This  nonlinear  formulation  also  makes  it  impossible  to  precompute  the 
Kalman  gains  and  covariance  matrices  (K  and  P),  a  technique  often  used  to  reduce 
filter  processing  time  for  linear  Kalman  filters.  Since  this  project  involves  post- 
flight  data  processing,  the  time  required  to  calculate  these  values  does  not  present 
a  major  problem. 

The  nonlinear  extended  Kalman  filter  requires  the  matrix  F,  which  consists  of 
the  partial  derivatives  of  Equation  (26)  with  respect  to  each  state  variable  evaluated 
at  the  current  state. 
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°. 

the  three  terms  Nx ,  Ny  and  Nz  in  the  tenth  row  of  the  F  matrix  in  Equation  (59) 
are  given  by: 


JV, 

Ny 

Nz 


1  dh 

T  dPt\ 

p*.p„p. 

1  dh 

T  OPy 

Px,P„P. 

i  dh 

t  dPz 

P*,P„P. 

(60) 

(61) 

(62) 


Nx ,  Ny,  and  Nz  are  the  partial  derivatives  of  altitude  with  respect  to  Px,  Py ,  and 
Pz  evaluated  at  the  current  aircraft  position.  Calculating  F  with  these  partial 
derivatives  gives  the  form  for  an  extended  Kalman  filter.  The  equation  for  aircraft 
altitude  given  position  is: 


Pz _ q(l  ~  e2) 

sin  L  Vl  —  e2  sin2  L 
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(63) 


This  equation  for  altitude  requires  the  aircraft  geodetic  latitude,  L.  The  geodetic 
latitude,  £«.,  is  not  equal  to  the  geocentric  latitude.  The  geodetic  latitude  is  given 
by  the  trigonometric  relationship 

pz 

Le  =  arcsin  . --.f - - -  (64) 

y/PFTPfTJ?  v  ' 

Geodetic  latitude  L  is  related  to  geocentric  latitude  Lc  using  the  relationship 


L  =  Le  +  f  sin(2£)  (65) 

where  /  is  the  flattening  ratio,  or  ellipticity,  based  on  the  semi-major  and  semi- 
minor  axes  of  the  ellipsoid  (4).  This  equation  cannot  be  solved  for  L  in  closed 
form  given  Le ,  but  it  can  be  expanded  into  a  sequence,  successively  evaluating  the 
right  side  of  the  equation  as  a  new  estimate  for  L.  When  treating  the  equation 
for  L  iteratively,  the  value  quickly  converges  to  within  0.001  degrees  of  the  proper 
answer.  Unfortunately,  convergence  to  the  exact  value  is  slow.  An  approximation 
is  made  by  including  three  terms  of  this  sequence. 

L  «  £c  +  / sin{2[£c  +  /  sin(2Ic)]}  (66) 

Knowing  the  relationship  between  L  and  the  position  vector  allows  computation 
of  the  partial  derivatives  with  respect  to  geodetic  latitude: 


—  e2)e2  sin  L  cos  L 

dL  ) 

(1  -  e2  sin2  L) 

1 

2 

dPt  ) 

The  partial  derivatives  of  L  can  be  taken  with  respect  to  the  geocentric 
latitude  using  the  relationship  between  Le  and  L  derived  in  Equation  (66): 


8L 

0Pt 


dL 
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dL 
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dL  dLc 
dLe  8PX 

[£.  +  /  »in{2(L,  +  f  sin(2Le)]}]  ~ 

(1  4-  2/cos{2£c  +  2/sin(2£c)} 

+  4/2  cos(2Ze)cos{2Z,e  +  2/sin(2Xc)}] 

[1  +  2/  cos{2£c  +  2/  sin(2Z,c)} 

dL 

+  Af2  cos(2 Le)  cos{2 Le  -f  2/  sin(2Ic)}] 
[1+2/  cos{2£c  +  2/  sin(2£c)} 

dL 

+  4/2  cos(2Le)cos{2Le  +  2/ sin(2£c)}]  -~ 


(70) 


(71) 


(72) 


Finally,  the  partial  derivatives  of  geocentric  latitude,  Le,  with  respect  to  Px, 
Py,  and  Pz  are  given  by: 


dLc 

dPx 


dLc 

dPv 

dLc 

dPz 
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dPx 


arcsm 


P ' 


.y/P:  +  PS+'P})\ 


-PxPz 

(p2  +  pv2  +  p^p2  +  p2 

_ -PyPz _ 

(P>  +  P*  +  P*)y/P*  +  PI 

Jpi^pl 

PI  +  P2  +  PI 


(73) 

(74) 

(75) 


Substituting  Equations  (73)  through  (75)  into  the  partial  derivatives  of  geode¬ 
tic  latitude,  Equations  (70)  through  (72),  and  substituting  those  results  into  the 
equations  for  Nx,  Ny  and  Nz  gives  the  required  partial  derivatives  to  fill  the  F 
matrix. 


Since  the  measurement  model  involves  nonlinear  functions  of  state  variables, 
the  filter  again  uses  partial  derivatives  to  calculate  matrix  elements  for  H.  The 
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required  range  partials  correspond  to  unit  line  of  site  vectors  and  are  given  by: 


Qj  1  CO 

1  II 

d 

dx 

sJ(P*  -Xb)2  +  (Pv  -  Yb)2  -f-  ( Pt  -  Zb)2 
P.-Xb 

(76) 

v/(ft  -  XbY  +  ( ft  -  1b)1  +  (ft  -  Zb)1 

dR 

Py-Ys 

(77) 

dy 

J(P r  -  XB)2  +  (Pv  -  Yb)2  +  ( P ,  -  ZB)2 

dR 

P'-Zb 

(78) 

dz 

,J(P*  -Xb)2  +  (Py  -  Yb)2  +  ( P ’  -  ZB)2 

Knowing  the  north  and  east  components  of  the  vector  from  the  beacon  to 

the  aircraft,  represented  by  N  and  E ,  along  with  the  magnetic  variation  (var)  at 
the  beacon  site,  allows  calculation  of  the  partial  derivatives  of  bearing  with  respect 
to  aircraft  position: 
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(79) 

(80) 
(81) 


Substituting  the  partial  derivatives  given  below  into  Equations  (79)  through  (81) 
allows  for  complete  evaluation  of  the  terms  required  in  the  filter: 


dE 

dx 

dE_ 

dy 

dE 

dz 

dN 

dx 


~{-sinlB(Px  ~  XB)  +  cos lB(Py  -  Yb)} 
dx 

—  sin  lB 
cos  lB 
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Q 

—  {-sin  LBcoslg(Px  -  XB ) 
dx 
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dN 

dy 

dN 

dz 


—  sin  Lb  sin  Ib{Pv  -  Vs)  4-  cos  Lb{Pz  -  Zb)} 

—  sin  Lb  cos  Ib 

—  sin  Lb  sin  4 

cos  Lb 


(85) 

(86) 

(87) 


4.2  System  Propagation 

The  operation  of  a  Kalman  filter  normally  separates  into  two  processes.  Be¬ 
tween  incorporating  measurements,  the  filter  propagates  the  last  estimate  of  the 
system  state  toward  the  next  measurement  time.  It  propagates  forward  in  time  for 
a  standard  filter,  and  both  forward  and  backward  in  time  for  an  optimal  smoother, 
which  is  covered  later  in  this  chapter.  Since  the  system  model  considered  here  has 
no  commanded  inputs,  the  term  u(f)  becomes  zero  in  the  general  filter  equations. 

Although  the  matrix  F  varies  with  time,  the  changes  occur  slowly.  This 
behavior  was  approximated  by  treating  F  as  piecewise  constant  over  the  0.5  sec 
sample  periods.  The  state  transition  matrix  to  propagate  the  system  states  forward 
in  time,  *£,  is  given  by: 


=  z:-1  {[si- F]-1} 


(88) 


0 

-1 

0 

0 

0 

0 

0 

0 

00000 

0 

0 

-1 

0 

0 

0 

0 

0 

00000 

3 

0 

0 

-1 

0 

0 

0 

0 

00000 

0 

s 

0 

0 

-1 

0 

0 

0 

00000 

0 

0 

3 

0 

0 

-1 

0 

0 

00000 

0 

0 

0 

3 

0 

0 

-1 

0 

00000 

0 

0 

0 

0 

3  —  r 

0 

0 

0 

00000 

0 

0 

0 

0 

0 

3  —  T 

0 

0 

00000 

0 

0 

0 

0 

0 

0 

s  —  r 

0 

00000 

■Nt 

0 

0 

0 

0 

0 

0 

3  — 

a  0  0  0  0  0 

0 

0 

0 

0 

0 

0 

0 

0 

5  0000 

0 

0 

0 

0 

0 

0 

0 

0 

05  000 

0 

0 

0 

0 

0 

0 

0 

0 

00500 

0 

0 

0 

0 

0 

0 

0 

0 

0005  0 

0 

0 

0 

0 

0 

0 

0 

0 

0000  5 

The  computer  program  MACSYMA  [54]  provided  the  inverse  of  this  matrix: 
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Table  3  lists  the  denominators  found  in  the  Laplace  transform  of  the  state 


transition  matrix.  These  have  been  listed  separately  due  to  the  complexity  of  the 
matrix. 


Table  3.  Matrix  Denominator  Terms 
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Table  4.  Laplace  Transform  Pairs 
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Table  4  gives  the  required  inverse  Laplace  transforms  for  this  state  transition 
matrix. 

4-3  System  Updates 

After  propagating  a  best  estimate  of  the  current  state  vector,  x(t,7 ),  and  the 
covariance  describing  the  uncertainty  of  that  value,  P {t~ ),  the  filter  incorporates 
available  measurements  to  improve  the  accuracy  of  those  estimated  values.  The 
Kalman  gain,  K(ij),  tells  the  filter  how  much  confidence  to  place  in  the  new  mea¬ 
surements,  relative  to  the  output  of  the  propagation  cycle.  For  extreme  Kalman 
gains,  the  filter  may  totally  ignore  its  own  propagated  state  estimate  and  place  full 
confidence  in  the  measurements,  or  the  filter  can  completely  ignore  the  measure¬ 
ments  if  it  does  not  consider  the  new  data  accurate  enough.  The  equations  used 
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to  perform  a  measurement  update  with  an  extended  Kalman  filter  are  given  by 
Equations  (54)  to  (56). 

4-4  Optimal  Smoothing 

The  Kalman  filter  equations  outlined  above  provide  the  optimal  estimate  of 
the  system  state  based  upon  all  measurements  up  to  and  including  the  one  at  time 
t{.  An  optimal  smoother  extends  that  idea  to  include  all  available  measurements, 
including  measurements  taken  in  the  future.  This  non-causal  nature  limits  optimal 
smoothing  to  post-flight  processing,  but  by  providing  additional  information  to  the 
estimator,  it  can  usually  provide  better  position  estimates  than  a  filter. 

4-4-1  Backward  Filter  The  optimal  smoother  uses  two  separate  Kalman  fil¬ 
ters.  The  forward  running  filter  uses  the  equations  from  the  preceding  section.  The 
second  filter  starts  at  the  final  time  and  runs  backward  toward  the  initial  time.  By 
combining  information  in  the  forward  filter,  incorporating  all  measurements  from 
initial  to  ti,  and  the  information  in  the  backward  filter,  which  uses  measurements 
U+i  to  tfinai,  the  smoother  can  use  all  the  measurements  to  estimate  the  system 
state. 

The  backward  filter  normally  runs  in  an  inverse  covariance  form  [35:6].  In 
this  form,  it  is  easier  to  work  with  a  state  vector  yj,  that  represents  the  product 
of  the  backward  inverse  covariance,  l,  and  the  backward  filter  state,  Xt,.  These 


values  are  initialized  with: 

y«.(*Hn.i)  =  «  (91) 

P*1^)  =  0  (92) 

For  each  measurement,  taken  at  time  ti,  the  backward  filter  is  updated  using: 

M‘t)  =  yi(<-)  +  Hr(i,)R-'(i,  )•((,)  (93) 

p;'((,+  )  =  P,-'((-)  +  Hr((i)R-'(<i)H(<i)  (94) 
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In  the  backward  filter,  t~  represents  the  instant  before  the  back  filter  receives 
the  measurement,  just  as  in  the  forward  filter.  However,  in  the  backward  filter,  t~ 
is  closer  to  ta„*\  than  and  tf  represents  the  instant  after  the.  measurement  is 

incorporated,  which  is  further  from  t final  than  from  < initial- 

Assuming  no  commanded  inputs,  the  backward  filter  propagates  to  the  next 
time  using: 

J(*)  =  Pr1(^)Gd(<i-i)[Gj(<,.1)P;1(tt)GJ(<i.1)  +  Q;H^i)]'1  (95) 


L(U)  =  I  —  J(tj)Gj(tj_i)  (96) 

M  (U)  =  L(U)P;l(tt)LT(U)  +  (97) 

ftftli)  =  *r(*,«i-»)L(*i)ft(t<+)  (98) 

=  ^T(t,>  t,_i)M(ti)$(tj,ti_i)  (99) 


In  the  preceding  equations,  G<j  is  the  identity  matrix,  and  Qj  corresponds 
to  the  integral  term  in  Equation  50: 

CM*,-i)=  P  -t[ti,<i-1;x(r)]G(r)Q(r)Gr(r)$r[<i,ti_1;x(r)](iT  (100) 

Jti-i 

4-4-2  Combining  Filter  Data  The  backward  filter  executes  after  the  forward 
filter  generates  state  and  covariance  estimates  for  the  entire  flight.  The  forward 
filter  uses  an  extended  Kalman  filter  with  a  nonlinear  system  model.  This  allows 
the  backward  filter  to  operate  with  a  linear  perturbation  filter  which  uses  the 
forward  filter  results  as  the  nominal  state  values.  For  each  measurement,  the 
smoother  calculates  state  estimates  by  combining  the  x("t* )  and  P (tf )  values  from 
the  forward  filter  with  the  y (,(/~)  and  P^t,-)  from  the  backward  filter  via: 


X(u) 

=  [i+piwonp 

(101) 

W(f.) 

=  p(!nxr(/i) 

(102) 

Y  (U) 

P(t«/tfin«a) 

=  I- W(i,)P,-'((-) 

=  Y(U)P(t*)YT(ti) 

(103) 

(104) 

(105) 


+w(<i)pl-,(<r)wr(ii) 

x(li/(«i»i)  =  X(l1)x((l')  +  P((i/tftaj)yk(i") 

4.5  Programming  Environment 

The  programs  written  for  this  project  are  coded  in  Fortran-77  on  a  Digital 
Equipment  Corporation  (DEC)  computer  running  the  VMS  operating  system.  The 
factors  supporting  this  programming  environment  are: 

•  John  Franzen’s  initial  programs  investigating  this  type  of  tracking  for  the 
4950th  Test  Wing  are  written  in  Fortran. 

•  The  data  analysts  in  the  4950th  Test  Wing  normally  use  Fortran  for  all  pro¬ 
gramming.  Using  it  for  this  project  makes  program  updates  and  maintenance 
easier. 

•  Both  AFIT  and  the  Test  Wing  computers  support  this  configuration,  easing 
the  transition  from  development  at  AFIT  to  use  in  the  Wing. 

4-6  System  Modeling 

As  explained  in  Chapter  III,  the  model  used  in  the  Kalman  smoother  con¬ 
tains  fifteen  states  in  order  to  model  positions,  velocities  and  accelerations  in  three 
dimensions,  the  lag  in  the  altimeter,  and  biases  in  the  five  measurements.  Driving 
this  system  are  eight  independent  noise  sources  which  enter  the  system’s  accelera¬ 
tion  and  bias  states  directly. 

4-6.1  Filter  Tuning  In  the  filter’s  model  of  the  real  world,  values  are  as¬ 
signed  to  the  strengths  of  the  driving  noise,  Q,  and  the  covariance  of  the  mea¬ 
surement  noise,  R.  Table  5  shows  the  values  chosen  for  the  diagonal  elements 
of  the  matrix  Q.  The  acceleration  driving  noise  strengths  correspond  to  the  3 <r 
value  for  aircraft  acceleration  of  15m/s2,  or  approximately  1.5y's  in  any  direction. 
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Table  5.  Driving  Noise  Strengths 


Ell 
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Accelerations 

ES9 

CM 

10-4m2/s 

10-flradJ/j 

TACAN  1  Range 
TAG  AN  1  Bearing 

Qw 

Qtj 

10"4m2/s 
10~erad3 / s 

TACAN  2  Range 
TACAN  2  Bearing 

100m3/ s 

Altimeter 

Table  6.  Measurement  Noise  Variances 


Cf 

Cal 

BiH 

Range  Bias  1 
Bearing  Bias  1 

ESQ 

E39 

400m2 

10“4rad2 

Range  Bias  2 
Bearing  Bias  2 

EH 

■m 

Altimeter  Bias 

These  values  for  the  first  three  diagonal  elements  of  Q  is  based  on  the  relationship 
Q  =  crT / 2  [37],  using  the  time  constant  T  from  Equation  (3).  Small  range  and  al¬ 
timeter  bias  driving  noise  strengths  reflect  the  essentially  constant  nature  of  these 
biases  while  preventing  the  tracker  from  totally  ignoring  future  measurements  re¬ 
lated  to  these  states.  These  values  are  the  result  of  iteratively  tuning  the  filter 
until  it  could  track  a  maneuvering  target. 

Table  6  shows  the  values  assigned  to  the  diagonal  elements  of  the  R  matrix. 

These  values  are  based  on  statistical  data  collected  by  the  4950th  Test  Wing  [39] 
and  the  Royal  Aircraft  Establishment  [46]. 

In  both  the  Q  and  R  matrices,  the  off-diagonal  terms  are  set  to  zero  which 
represents  statistical  independence  of  the  measurements.  This  assumes  errors  from 
one  beacon  are  independent  of  errors  from  a  second  beacon  since  they  are  geograph¬ 
ically  isolated  from  each  other.  The  range  and  bearing  information  provided  by  ( 
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a  single  TACAN  are  considered  independent  since  they  involve  different  measure¬ 
ment  mechanisms.  This  same  reasoning  allows  treating  the  barometrically  mea¬ 
sured  altitude  as  independent  of  the  ranges  and  bearings  determined  with  radio 
signals. 
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V.  Computer  Simulations 


A  complete  evaluation  of  the  tracking  technique  developed  in  this  thesis  re¬ 
quires  TACAN  and  altimeter  measurements  collected  while  a  high-precision  ref¬ 
erence  system  tracks  the  aircraft.  Funding  and  scheduling  considerations  have 
prevented  the  Test  Wing  from  providing  this  type  of  data.  Simulated  measure¬ 
ment  data  based  on  estimates  of  real-world  noise  characteristics  provide  the  basis 
for  this  chapter. 

5.1  Measurement  Data  Generation 

Input  data  to  the  dual  TACAN  tracker  is  generated  using  a  combination  of 
five  computer  programs.  Aircraft  maneuver  commands  are  entered  into  PROF- 
GEN,  an  aircraft  flight  profile  generating  program  developed  by  Stan  Musick  at 
the  Air  Force  Wright  Aeronautical  Laboratories  [40].  This  program  uses  a  high- 
order  numerical  integrator  to  solve  the  aircraft  equations  of  motion  and  provides 
realistic  values  for  aircraft  position,  velocity  and  acceleration. 

The  output  of  PROFGEN  is  combined  with  a  list  of  TACAN  beacons  used 
during  the  flight  and  the  TACAN  positions.  The  resulting  file  contains  the  true 
ranges  and  bearings  from  the  beacons  to  the  aircraft.  Finally,  the  measurement 
file  is  corrupted  by  noise,  with  each  TACAN  beacon  having  its  own  noise  statistics 
assigned  to  representative  values  (39,46].  Table  7  outlines  the  steps  necessary  to 
create  a  simulated  measurement  file. 

5.2  Test  Scenario 

The  measurement  file  simulates  a  15-minute  test  flight  starting  from  a  sta¬ 
tionary  position  on  the  runway  at  Wright-Patterson  AFB.  At  time  zero,  the  plane 
accelerates  and  takes  off  to  the  southwest.  After  climbing  to  an  altitude  of  1600 
meters,  the  plane  executes  a  slow  turn  to  the  left.  As  the  flight  progresses,  the 
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Table  7.  Creating  Simulation  Data 


I 


Input  Data 

Program 

Output  Data 

Maneuver  Commands 

PROFGEN 

Aircraft  Position 

Aircraft  Position 
TACAN  Positions 

INSTRM 

Noise  Free  Altimeter 
and  TACAN  Measurements 

Noise  Free 
Measurements 

CORRUPT 

Noise  Corrupted 
Measurements 

Figure  7.  Test  Flight  Ground  Track 


plane  also  performs  a  slow  turn  to  the  right  and  four  high-speed  turns.  The  sample 
flight  includes  several  velocity  and  altitude  changes.  These  maneuvers  represent  the 
type  of  flying  usually  involved  in  flight  tests  conducted  by  the  4950th  Test  Wing. 
Figure  7  gives  a  map-like  view  of  this  flight  path,  with  the  aircraft  travelling  from 
the  upper  right  to  the  lower  left  of  the  plot. 

Figures  8  and  9  show  the  aircraft  altitude  and  velocity  plotted  as  functions 
of  time.  Both  graphs  begin  with  the  aircraft  stationary  on  the  runway  at  time 
t  =  —40  seconds.  The  aircraft  begins  the  take-off  acceleration  at  i  =  0  and  takes 
off  at  i  =  15  seconds. 
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Figure  10.  TACAN  Beacon  Locations 

An  important  operation  for  the  tracker  involves  adjusting  to  constantly  vary¬ 
ing  measurement  sources  as  different  pairs  of  TACAN  beacons  provide  measure¬ 
ments  during  a  flight.  This  test  scenario  uses  six  different  beacons,  retuning  the 
on-board  transceivers  at  times  t  =  —10,  t  =  45,  t  =  200,  and  t  —  500  seconds. 
Figure  10  shows  the  position  of  the  aircraft  track  in  relation  to  the  ground  beacons. 

The  second  reason  for  changing  beacons  throughout  the  flight  involves  the 
uncertainty  involved  when  resolving  multiple  range  measurements.  The  minimum 
uncertainty  occurs  when  the  difference  in  bearings  from  the  aircraft  to  the  two 
TACAN  beacons,  or  cut  angle,  equals  90  degrees  [39:F3].  TACAN  selection  should 
maintain  this  cut  angle  between  60  and  120  degrees  for  best  tracker  performance. 
During  this  test  flight,  the  angle  between  the  TACANs  varies  from  40  to  85  degrees. 

5.3  System  Performance 

5.3.1  Single  Simulation  Performance  The  Kalman  smoother  performance 
is  evaluated  by  comparing  smoother  output,  expressed  in  earth-centered,  earth- 
fixed  coordinates  with  the  original,  noise-free  values  derived  from  PROFGEN. 
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Figure  11  shows  the  component  of  the  position  error  in  the  x  direction,  which 
is  calculated  by  subtracting  the  tracker  value  for  Pt  from  the  true  value  Px.  This 
graph  has  a  temporally  averaged  mean  error  of  —9  meters  and  a  temporally  derived 
standard  deviation  of  14  meters. 

This  same  technique  was  used  to  generate  the  error  plots  in  the  y  and  z 
directions  shown  in  Figures  12  and  13.  The  y  error  component  has  a  temporally 
averaged  mean  value  of  15  meters  and  a  temporal  standard  deviation  of  22  meters. 
The  z  component  of  error  has  a  temporally  averaged  mean  value  of  42  meters  and 
a  temporal  standard  deviation  of  21  meters. 

Radial  error  provides  another  means  to  assess  tracker  accuracy.  This  is  the 
method  of  evaluating  navigation  accuracy  outlined  in  the  North  Atlantic  Treaty 
Organization  (NATO)  Standardization  Agreement  (STAN AG)  4278,  “Method  of 
Expressing  Navigation  Accuracy”  [5j.  The  radial  error  distance  (Rf)  is  computed 
using  the  true  position  (xe,  yt,  zt )  and  the  smoother  estimate  of  position  ( Px ,  Py, 

A).  ^ _ 

Re  =  \!(*t  -  Px)2  +  (yt  -  Py)2  +  {zt  -  Pzf  ( 106) 
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_ TIME _ 

Figure  14.  Example  of  Radial  Error 

Figure  14  shows  the  radial  error  for  one  simulation  run  plotted  versus  time. 

The  shift  in  the  error  plot  at  time  t  =  200  seconds  corresponds  to  a  change  in 
TACAN  beacons  used  by  the  system.  The  opportunity  to  recalculate  measurement 
biases  and  more  favorable  solution  geometry  allows  the  filter  to  estimate  the  aircraft 
position  more  accurately. 

According  to  STAN  AG  4278,  a  positioning  system  should  have  a  radial  error 
less  than  the  quoted  radial  accuracy  95  percent  of  the  time.  Figure  15  plots  the 
cumulative  probability  for  radial  error  based  upon  the  simulation  data.  Following 
this  standard,  radial  accuracy  for  this  system  is  90  meters. 

The  complex  calculations  associated  with  a  Kalman  smoother  make  it  im¬ 
portant  to  compare  the  performance  of  the  smoother  against  a  single  pass  Kalman 
filter.  A  small  improvement  in  tracking  performance  may  not  justify  the  additional 
computational  loading  needed  to  smooth  the  data.  The  output  of  the  extended 
Kalman  filter  operating  forward  in  time  is  compared  with  the  final  smoothed  val¬ 
ues.  Figure  16  shows  the  cumulative  radial  error  distributions  for  both  the  filter 
and  the  smoother.  At  the  95  percent  level,  the  forward  filter  has  a  STANAG  4278  I 


43 


I 


I 


44 


Figure  17.  Monte  Carlo  Results  in  X  Direction 
accuracy  of  109  meters  compared  with  the  90  meter  accuracy  for  the  full  smoother. 

5.3.2  Monte  Carlo  Results  A  total  of  11  Monte  Carlo  runs  were  conducted. 
Each  Monte  Carlo  run  processes  a  different  corrupted  version  of  the  same  aircraft 
trajectory.  The  mean  and  standard  deviation  of  PT ,  Pv  and  Pz  are  calculated  for  the 
11  values  at  each  sample  period.  Figures  17-19  show  the  mean  plus  or  minus  one 
standard  deviation  throughout  the  flight.  The  jumps  that  occur  in  these  plots  at 
times  t  —  —10,  t  =  45,  t  =  200  and  t  =  500  seconds  correspond  to  TACAN  beacon 
changes.  Changing  TACAN  beacons  also  causes  changes  in  solution  geometry. 

Using  the  eleven  Monte  Carlo  values  for  radial  error  at  each  sample  period, 
the  mean  and  the  mean  plus  or  minus  one  standard  deviation  are  calculated  and 
plotted  in  Figure  20.  Figure  21  shows  the  cumulative  distribution  of  the  radial 
error  for  the  Monte  Carlo  runs.  The  95  percent  accuracy  figure  for  this  system  is 
135  meters. 

Changing  aircraft  position  along  with  changing  TACAN  beacon  geometry  was 
shown  earlier  to  cause  large  changes  in  the  smoother  performance.  This  bimodal 
nature  of  the  radial  error  cumulative  distribution  also  demonstrates  this  effect. 
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Considering  only  the  values  after  the  TACAN  change  at  time  t  —  200  seconds 
gives  the  cumulative  distribution  shown  in  Figure  22.  The  system  achieves  an 
accuracy  of  109  meters  from  t  =  200  seconds  through  the  end  of  the  flight  at 
t  ~  1000  seconds. 
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VI.  Conclusions  and  Recommendations 


6.1  System  Accuracy 

Performance  of  this  tracking  system  meets  the  desired  accuracy  specifica¬ 
tions.  During  the  computer  simulations,  the  system  placed  the  aircraft  within 
135  meters  of  its  true  position  95  percent  of  the  time.  The  system  does  not  have 
any  difficulty  tracking  the  aircraft  through  accelerations,  decelerations,  ascents, 
descents  or  turns. 

6.2  Recommendations 

The  tracking  system  developed  in  this  thesis  cannot  satisfy  the  requirements 
of  all  test  programs,  but  it  does  offer  a  low-cost  alternative  to  flying  at  an  instru¬ 
mented  test  range.  This  method  is  adequate  for  many  projects.  Other  programs 
requiring  more  accurate  tracking  information  can  use  this  system  to  supplement 
test  range  flying.  It  would  allow  thorough  equipment  and  procedural  checkouts  for 
complex  programs  before  spending  the  time  and  money  to  deploy  to  a  range  for 
full-scale  testing.  The  Test  Wing  should  follow  up  on  this  project  and  turn  the 
promise  shown  in  these  results  into  an  operational  system. 

6.3  Follow-On  Work 

After  working  on  this  project,  several  possible  extensions  have  surfaced  that 
will  improve  program  execution  and  enhance  systems  tracking  capability. 

6.3.1  Binary  Data  Files  The  program  currently  does  all  file  manipulations 
using  ASCII  data  files.  These  files  allow  easier  program  debugging  since  the  VAX 
text  editor  can  examine  or  modify  the  information. 

All  values  are  double  precision  and  take  up  20  bytes  of  storage  per  value. 
Using  binary  files  in  the  Fortran  code  to  read  and  write  to  disk  would  decrease 
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storage  requirements  to  8  bytes  per  value,  a  savings  of  60  percent.  Along  with 
storage  savings,  these  files  allow  faster  disk  reads  and  writes  since  fewer  bytes 
are  transferred,  and  the  computer  does  not  have  to  convert  between  ASCII  and 
internal  floating  representations. 

6.8.2  Adaptive  Extended  Kalman  Filter  The  current  tracking  system  uses  a 
15-state  extended  Kalman  filter.  Five  of  those  states  estimate  essentially  constant 
values.  A  10-state  adaptive  Kalman  filter  using  those  bias  values  as  parameters 
can  provide  similar  tracking  performance.  The  lower  dimensionality  of  the  ex¬ 
tended  filter  decreases  system  storage  requirements  and  allows  faster  processing  of 
measurements. 

6.3.3  Factored  Filter  The  current  filter  performs  full-state  estimation  using 
the  standard  equations  for  an  extended  Kalman  filter.  To  avoid  numerical  diffi¬ 
culties  that  can  arise  in  the  filter  equations,  all  calculations  use  double  precision 
variables  on  the  computer.  By  recoding  the  filter  using  the  UD  factored  form,  most 
of  these  operations  can  be  done  using  single  precision  math  [34:392]. 

6.3.4  Fine  Tuning  The  lack  of  dual  TACAN  data  collected  in  the  real  world 
does  not  make  it  worthwhile  to  expend  a  lot  of  effort  in  tuning  the  Q  and  R  matrices 
used  by  the  Kalman  filter.  When  the  Test  Wing  has  this  data  available,  then  the 
filter  should  be  tuned  to  make  allowances  for  unmodelled  effects  that  may  corrupt 
the  measurement  signals.  Another  approach  to  tuning  the  filter  would  involve 
developing  higher  order  truth  models  that  better  represent  real  world  data. 

6.3.5  Integration  With  Winy  Data  Collection  and  Reduction  Test  data  for¬ 
mats  often  vary  from  one  program  to  the  next.  Working  with  the  Test  Wing  data 
analysts  and  instrumentation  personnel  to  define  standards  for  using  this  system 
will  avoid  wasted  time  and  effort  in  the  future.  These  standards  should  include: 
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•  an  equipment  package  for  installation  on  the  test  aircraft 

•  guidelines  for  selecting  TACAN  beacons  during  a  test  mission 

•  in-flight  operating  procedures  for  the  TACAN  transceivers 

•  automating  the  procedures  for  converting  raw  dual  TACAN  data  into  track¬ 
ing  data 


6.3.6  Multiple  Sensor  Integration  This  research  effort  concentrated  on  the 
measurements  available  in  a  minimally  equipped  aircraft.  Many  Wing  test  pro¬ 
grams  use  additional  on-board  navigation  equipment,  such  as  INS,  GPS  and  long 
range  navigation  (LORAN).  Developing  a  Kalman  filter  or  filter  to  integrate  these 
various  measurement  sources  will  allow  development  of  a  more  complex,  and  hope¬ 
fully  more  accurate,  tracker. 
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